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Abstract 

Background: Herbivore grazing is a multiple-component process that includes wounding, defoliation, and saliva deposition. 
Despite the extensive published research on mechanical wounding and defoliation, no analysis to identify the genes that 
specify defoliation and mechanical wounding has been performed. Moreover, the influence of the expression of these 
genes on plant regrowth after defoliation remains poorly understood. 

Results: Seven cDNA libraries for RNA samples collected from stubble tissues that had been mechanically wounded or 
defoliated at 2, 6 and 24 h along with the control were sequenced using the Illumina/Solexa platform. A comparative 
transcriptomic analysis of the sequencing data was conducted. In total, 1,836 and 3,238 genes were detected with 
significant differential expression levels after wounding and defoliation, respectively, during one day. GO, KOG and 
pathway-based enrichment analyses were performed to determine and further understand the biological functions of those 
differentially expressed genes (DEGs). The results demonstrated that both wounding and defoliation activated the systemic 
synthesis of jasmonate (JA). However, defoliation specifically reduced the expression levels of ribosomal protein genes, cell 
division or cell expansion-related genes, and lignin biosynthesis genes and may have negatively affected plant growth. 
Further analysis revealed that the regrowth of elongating leaves was significantly retarded after defoliation at 6 h through 
the following 7 days of measurement, suggesting that the gene expression pattern and phenotype are consistent. Fifteen 
genes were selected, and their expression levels were confirmed by quantitative RT-PCR (qRT-PCR). Thirteen of them 
exhibited expression patterns consistent with the digital gene expression (DGE) data. 

Conclusions: These sequencing datasets allowed us to elucidate the common and distinct mechanisms of plant responses 
to defoliation and wounding. Additionally, the distinct DEGs represent a valuable resource for novel gene discovery that 
may improve plant resistance to defoliation from various processes. 
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Introduction 

Grasslands are a major part of the global ecosystem, covering 
37% of the Earth's terrestrial area [1]. Herbivory is a fundamental 
driver of grassland plant community composition. Nearly all 
grasslands worldwide have been grazed by populations of wild and 
domestic herbivores, such as cattle [2], horses [3], sheep [4], deer 
[5], rabbits [6], and kangaroos [7]. Herbivore grazing is a 
multiple-component process that includes wounding, defoliation, 
and saliva deposition [8]. 

Wounding responses have been demonstrated in species 
throughout the plant kingdom, and previous research has focused 



on identifying systemic wound signals and the mechanisms by 
which they are generated, transported, and perceived [9]. Leaf 
wounding is a well-characterized stimulus for the rapid activation 
of jasmonate (JA) synthesis through distinct autonomous and non- 
autonomous cellular pathways [9-12]. 

Defoliation not only induces serious wounding but also removes 
a large fraction of photosynthetically active leaves from grass 
plants. Therefore, plant photosynthesis is reduced; carbohydrate 
reserves stored in the plant are mobilized; and the transcript levels 
of photosynthesis-related and carbohydrate metabolism-related 
genes, such as genes encoding the small subunits of Rubisco, 
Rubisco activase, light harvesting chlorophyll a/b binding 
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Figure 1. Summary of differentially expressed genes. The categories of DEGs expressed at one time point and at overlapping time points for 
wounding and defoliation at 2, 6, and 24 h. (A) A summary of the up- and down-regulated DEG numbers. (B) and (C) DEGs expressed at one time 
point and overlapping time points for wounding and defoliation, respectively. 
doi:1 0.1 371 /journal.pone.0089495.g001 



proteins, fructan exohydrolase (FEH), sucrose transporter (SUT), 
invertase, P-amylase, and starch synthase, are changed [13-18]. 

We previously identified 466 genes that responded to simulated 
grazing using the Affymetrix Rice GeneChips. These genes encode 
proteins that participate in signal transduction, miRNA regulation, 
cell wall modification, carbohydrate metabolism, hormone 
synthesis, and molecular transport [8], However, it is unclear 
whether the responses were due to wounding, defoliation, or other 
factors. 



Currently, the genes that specify defoliation and mechanical 
wounding remain unknown, and little is understood about how 
expression of these genes influences plant regrowth after defoli- 
ation. In this study, the transcriptomes of sheepgrass stubble tissues 
collected from defoliated and mechanically wounded plants were 
compared using the Solexa sequencing system. By comparing the 
differentially expressed gene (DEG) responses to defoliation and 
wounding, the common and distinct pathways and genes between 
defoliation and wounding were identified. These sequencing 
datasets allowed us to elucidate the common and distinct 
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Figure 2. Histogram of functional Gene Ontology analysis of DEGs. The frequency of GO terms was analyzed using the GO Slim Assignment. 
The x- and y-axes represent the names of clusters and the ratio of each cluster, respectively. W, mechanical wounding; D, defoliation. 
doi:1 0.1 371 /journal.pone.0089495.g002 



mechanisms of the plant response to defoliation and wounding, 
and the distinct DEGs represent a valuable resource for novel gene 
discovery that may improve plant tolerance toward defoliation 
caused by herbivores, insects, pathogens, or mechanical stress. 

Materials and Methods 

Plant Materials, Growing Conditions, and Treatments 

The sheepgrass variety Zhongke 3 was grown in a soil mix of 
peat moss and vermiculite (2:1, v/v) in a greenhouse at 23°C with 
a photoperiod of 16 h light/8 h dark. Wounding and defoliation 
treatments were performed on eight-week-old sheepgrass seed- 
lings. To induce defoliation, approximately two-thirds of the 
aboveground biomass was removed, and to induce wounding, 
tweezers were used to mechanically wound a site at the same 
position as the defoliation site. The remaining one-third of the 
aboveground biomass was collected 2, 6, and 24 h after 
defoliation, and the corresponding parts were collected for control 
seedlings and wounding at 2, 6, and 24 h. Each sample included 3 
replicate pots and 21 seedlings. The samples were collected and 
immediately frozen in liquid nitrogen and stored at — 80°C. 

RNA Extraction and Solexa/lllumina Sequencing 

Total RNA was extracted from individual control, wounded 
and defoliated plants at 2, 6, and 24 h using TRIzol reagent 
(Invitrogen, Carlsbad, CA, USA). After precipitation, the RNA 
was purified using Qiagen's RNeasy Kit with on-column DNase I 
digestion according to the manufacturer's instructions. RNA 
concentrations were measured using the Bioanalyzer 2 1 00 (Agilent 
Technologies, Inc., Waldbronn, Germany), and the RNA integrity 
was analyzed on a 1.0% (w/v) agarose gel. Approximately 5 to 
8 (Xg total RNA was used to construct each RNA-seq library. Poly 
(A) mRNAs were isolated from the total RNA using oligo (dT) 



magnetic beads (Illumina, San Diego, CA). RNA fragmentation, 
cDNA synthesis, and PCR amplification were performed accord- 
ing to the Illumina RNA-seq protocol (Cat # RS-100-0801). 
Seven cDNA libraries, which were isolated from mechanically 
wounded and defoliated material at 2, 6, and 24 h and the control, 
were sequenced using the Illumina HiSeq 2000 System at the 
Chinese National Human Genome Center (Shanghai, China). 
The sequence reads generated in this study have been deposited in 
the NCBI sequence read archive (SRA065691). 

Analysis and Screening of Differentially Expressed Genes 
(DEGs) 

The reads for each sample were counted and mapped back to 
previously generated sheepgrass reference genes [19] using the 
software tool Bowtie [20]. The read number of each gene was 
transformed into RPKM (Reads Per Kilobases per Million reads) 
[21], and the differendy expressed genes were identified by a 
DEGseq package using the MARS method (MA-plot-based 
method with Random Sampling model) [22]. The p-value was 
obtained from the differential gene expression test. FDR (False 
Discovery Rate) manipulation was used to determine the p-value 
threshold in multiple tests and analyses. Both FDR <0.001 and 
the absolute value of | the log2Ratio Si were used as thresholds 
to identify the significant DEGs [23]. 

Gene Ontology, EuKaryotic Orthologous Groups and 
Pathway Enrichment Analysis 

Gene Ontology (GO) is an international standardized functional 
gene classification system that describes the properties of genes and 
their products in any organism using three ontologies: cellular 
component, molecular function, and biological process. GO 
enrichment analysis was performed using agriGO (http:// 
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Table 1. KOG enrichment analysis of the wounding-responsive DEGs compared with the reference transcriptome. 



KOG functional cluster 


Transcriptome frequency of use 


DEGs frequency of use 


p-value 


g-value 


Chromatin structure and dynamics 


636 of 32382 (1.96%) 


85 of 832 (10.22%) 


7.16E-37 


1.79E-35 


Energy production and conversion 


1 062 of 32382 (3.28%) 


46 of 832 (5.53%) 


0.0002 


0.003 


Amino acid transport and metabolism 


1 385 of 32382 (4.28%) 


54 of 832 (6.49%) 


0.001 


0.0085 


Carbohydrate transport and metabolism 


1 897 of 32382 (5.86%) 


62 of 832 (7.45%) 


0.02 


0.14 


Defense mechanisms 


333 of 32382 (1.03%) 


12 of 832 (1.44%) 


0.09 


0.43 


Secondary metabolites biosynthesis, transport and catabolism 


1511 of 32382 (4.67%) 


46 of 832 (5.53%) 


0.1 


0.43 


Nucleotide transport and metabolism 


398 of 32382 (1.23%) 


13 of 832 (1.56%) 


0.15 


0.53 


Cell wall/membrane/envelope biogenesis 


424 of 32382 (1.31%) 


13 of 832 (1.56%) 


0.21 


0.64 


Inorganic ion transport and metabolism 


942 of 32382 (2.91%) 


27 of 832 (3.25%) 


0.24 


0.67 


Extracellular structures 


1 5 of 32382 (0.05%) 


0 


0.32 


0.77 


Cell motility 


1 6 of 32382 (0.05%) 


0 


0.34 


0.77 


Transcription 


1 520 of 32382 (4.69%) 


40 of 832 (4.81%) 


0.39 


0.83 


Cytoskeleton 


1020 of 32382 (3.15%) 


26 of 832 (3.13%) 


0.46 


0.88 


Lipid transport and metabolism 


1 228 of 32382 (3.79%) 


31 of 832 (3.73%) 


0.49 


0.88 


Translation, ribosomal structure and biogenesis 


1819 of 32382 (5.62%) 


44 of 832 (5.29%) 


0.63 




Posttranslational modification, protein turnover, chaperones 


3144 of 32382 (9.71%) 


77 of 832 (9.25%) 


0.65 


^^^^ 


Coenzyme transport and metabolism 


322 of 32382 (0.99%) 


4 of 832 (0.48%) 


0.92 




Nuclear structure 


142 of 32382 (0.44%) 


0 


0.98 




Replication, recombination and repair 


821 of 32382 (2.54%) 


12 of 832 (1.44%) 


0.98 




Intracellular trafficking, secretion, and vesicular transport 


1589 of 32382 (4.91%) 


28 of 832 (3.37%) 


0.98 




Signal transduction mechanisms 


51 14 of 32382 (15.79%) 


107 of 832 (12.86%) 


1 




Cell cycle control, cell division, chromosome partitioning 


683 of 32382 (2.11%) 


8 of 832 (0.96%) 


1 




General function prediction only 


3292 of 32382 (10.17%) 


63 of 832 (7.57%) 


1 




Function unknown 


1685 of 32382 (5.20%) 


21 of 832 (2.52%) 


1 




RNA processing and modification 


1 384 of 32382 (4.27%) 


13 of 832 (1.56%) 


1 
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bioinfo.cau.edu.cn/agriGO/, [24]). EuKaryotic Orthologous 
Groups (KOG) enrichment analysis was conducted through 
hypergeometric distribution testing using the Phyper function in 
the R software package (http://www.rproject.org/). The Bonfer- 
roni correction was used to adjust the p-values. The significandy 
enriched functional clusters were selected based on a corrected q- 
value less than 0.05. Kyoto Encyclopedia of Genes and Genomes 
(KEGG) pathway enrichment analysis was conducted using 
KOBAS 2.0 (http://kobas.cbi.pku.edu.cn/, [25]). 

Real-time Quantitative RT-PCR (qRT-PCR) Analysis 

Real-time quantitative RT-PCR (qRT-PCR) analysis was used 
to verify the DGE results. The RNA samples used for the qRT- 
PCR assays were the same as in the DGE experiments. Gene- 
specific primers were designed according to the reference unigene 
sequences using Primer Premier 6.0. Fifteen genes were selected 
from the DEGs under defoliation at 2 h for quantitative RT-PCR 
assays, and these genes were involved in JA biosynthesis, the 
gibberellin signaling pathway, cell wall biogenesis/degradation, 
and starch and sucrose metabolism (Table SI). qRT-PCR was 
performed according to the manufacturer's specifications (Ta- 
KaRa SYBR Premix Ex Taq II Kit, Dalian, China). The following 
SYBR Green PCR cycling conditions were used: denaturation at 
95°C for 10 s, followed by 40 cycles of 94°C for 5 s and 60°C for 
20 s. The PCR experiments were performed using an iO_ 5 
Multicolor Real-time PCR Detection System (Bio-Rad, USA). 



The results were normalized to the sheepgrass actin gene, and the 
relative gene expression levels were calculated using 2~ CT . 

Results 

Wounding and Defoliation Exhibit Distinct DEGs 

To identify the DEGs involved in wounding and defoliation, we 
sequenced seven cDNA libraries and generated 192 million clean 
reads ranged from 4,811,093 to 6,411,192 per library after 
filtering the low-quality and adapter sequences. These clean reads 
were mapped to the sheepgrass reference genes as described in the 
Materials and methods. In total, 36,413, 35,925, 35,113, 35,929, 
34,534, 36,828, and 36,956 unigenes were detected for sequenced 
libraries 1-7, respectively (Table S2). The reference unigenes 
mapped by all clean reads were enhanced as the sequencing 
amount (total read number) increased. However, when the 
sequencing counts reached 1.5 million reads or higher, the 
number of detected genes became saturated (Figure SI), indicating 
that the sequencing depth was sufficient to cover the transcrip- 
tome. 

By comparing the treated and control Solexa libraries, a great 
number of differentially expressed transcripts were identified. To 
study the subset of genes associated with wounding and 
defoliation, we analyzed the most differentially regulated genes 
with a | log2Ratio | > 1 and a FDR <0.001. There were 706 (356 
up-regulated and 350 down-regulated), 891 (442 up-regulated and 
449 down-regulated), and 803 (362 up-regulated and 441 down- 
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Table 2. KOG enrichment analysis of the defoliation-responsive DEGs compared with the reference transcriptome. 



KOG functional cluster 


Transcriptome frequency of use 


DEGs frequency of use 


p- value 


g-value 


Chromatin structure and dynamics 


636 of 32382 (1.96%) 


144 of 1605 (8.97%) 


9.95E-56 


2.49E-54 


Carbohydrate transport and metabolism 


1897 of 32382 (5.86%) 


156 of 1605 (9.72%) 


1.53E-10 


1.92E-09 


Cytoskeleton 


1020 of 32382 (3.15%) 


73 of 1605 (4.55%) 


0.0007 


0.0045 


Translation, ribosomal structure and biogenesis 


1819 of 32382 (5.62%) 


126 of 1605 (7.85%) 


5.97E-05 


0.0005 


Inorganic ion transport and metabolism 


942 of 32382 (2.91%) 


59 of 1605 (3.68%) 


0.03 


0.14 


Secondary metabolites biosynthesis, transport and catabolism 


1511 of 32382 (4.67%) 


87 of 1605 (5.42%) 


0.06 


0.27 


Cell wall/membrane/envelope biogenesis 


424 of 32382 (1.31%) 


27 of 1605 (1.68%) 


0.07 


0.27 


Defense mechanisms 


333 of 32382 (1.03%) 


21 of 1605 (1.31%) 


0.1 1 


0.29 


Nucleotide transport and metabolism 


398 of 32382 (1.23%) 


25 of 1605 (1.56%) 


0.09 


0.29 


Extracellular structures 


15 of 32382 (0.05%) 


1 of 1605 (0.06%) 


0.17 


0.42 


Energy production and conversion 


1 062 of 32382 (3.28%) 


58 of 1605 (3.61%) 


0.19 


0.45 


Transcription 


1 520 of 32382 (4.69%) 


77 of 1605 (4.80%) 


0.39 


0.78 


Nuclear structure 


142 of 32382 (0.44%) 


7 of 1605 (0.44%) 


0.41 


0.78 


Lipid transport and metabolism 


1228 of 32382 (3.79%) 


61 of 1605 (3.80%) 


0.46 


0.82 


Cell motility 


16 of 32382 (0.05%) 


o 


0.56 


0.93 


Posttranslational modification, protein turnover, chaperones 


3144 of 32382 (9.71%) 


152 of 1605 (9.47%) 


0.61 


0.95 


Amino acid transport and metabolism 


1 385 of 32382 (4.28%) 


56 of 1605 (3.49%) 


0.94 




Coenzyme transport and metabolism 


322 of 32382 (0.99%) 


12 of 1605 (0.75%) 


0.81 


^^^^ 


Replication, recombination and repair 


821 of 32382 (2.54%) 


23 Of 1605 (1.43%) 






Intracellular trafficking, secretion, and vesicular transport 


1589 of 32382 (4.91%) 


50 of 1605 (3.11%) 






Signal transduction mechanisms 


51 14 of 32382 (15.79%) 


159 of 1605 (9.91%) 






Cell cycle control, cell division, chromosome partitioning 


683 of 32382 (2.11%) 


18 of 1605 (1.12%) 






General function prediction only 


3292 of 32382 (10.17%) 


120 of 1605 (7.48%) 






Function unknown 


1685 of 32382 (5.20%) 


49 of 1605 (3.05%) 






RNA processing and modification 


1 384 of 32382 (4.27%) 


44 of 1605 (2.74%) 
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regulated) genes that responded to wounding at 2, 6, and 24 h, 
respectively. For defoliation, we identified 2,21 1 (685 up-regulated 
and 1526 down-regulated), 997 (451 up-regulated and 546 down- 
regulated), and 1,111 (480 up-regulated and 631 down-regulated) 
genes that were significantly expressed at 2, 6, and 24 h, 
respectively (Figure 1A). These significantly differentially wound- 
ing- and defoliation-responsive genes were further divided into 
genes expressed at one time point and those expressed at 



overlapping time points. For example, 1,836 and 3,238 genes 
were differentially expressed between the wounding and defolia- 
tion treatments, respectively, during one day (Figure IB, C), and 
954 genes were uniquely expressed under only defoliation. These 
DEGs are listed in Table S3 and Table S4, respectively. 



Table 3. KEGG pathways that are enriched in the wounding-responsive DEGs compared with the reference transcriptome. 



KEGG pathway 


DEGs number 


Background number 


P-Value 


Corrected P-Value 


Valine, leucine and isoleucine degradation 


19 


71 


5.28E-07 


3.24E-05 


Linoleic acid metabolism 


10 


23 


2.08E-06 


9.55E-05 


Ascorbate and aldarate metabolism 


14 


50 


9.35E-06 


0.000344057 


Starch and sucrose metabolism 


31 


188 


1 .67E-05 


0.000512626 


Fatty acid elongation 


8 


22 


0.0001044 


0.002745301 


Arginine and proline metabolism 


19 


102 


0.0001 393 


0.003203327 


alpha-Linolenic acid metabolism 


11 


46 


0.0003908 


0.00798909 


Lysine degradation 


9 


41 


0.0024838 


0.045702656 


Lysine biosynthesis 


6 


21 


0.0032349 


0.049957314 


Plant hormone signal transduction 


17 


112 


0.0032581 


0.049957314 



doi:1 0.1 371 /journal.pone.0089495.t003 
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Table 4. KEGG pathways that are enriched in the defoliation-responsive DEGs compared with the reference transcriptome. 



KEGG pathway 


DEGs number 


Background number 


/'-Value 


Corrected P-Value 


Starch and sucrose metabolism 


62 


188 


1.85E-10 


1 .06E-08 


Ribosome 


91 


323 


1.97E-10 


1 .06E-08 


alpha-Linolenic acid metabolism 


24 


46 


2.77E-09 


1.19E-07 


Phenylpropanoid biosynthesis 


28 


62 


8.87E-09 


3.18E-07 


Flavonoid biosynthesis 


15 


22 


1 .94E-08 


5.96E-07 


Stilbenoid, diarylheptanoid and gingerol biosynthesis 


9 


11 


1 .34E-06 


3.59E-05 


Linoleic acid metabolism 


13 


23 


3.97E-06 


9.48E-05 


Phenylalanine metabolism 


22 


62 


3.96E-05 


0.000850667 



doi:1 0.1 371 /journal.pone.0089495.t004 



GO Functional Analysis of DEGs 

A total of 835 DEGs that responded to wounding and 1,560 
DEGs that responded to defoliation were divided into 44 GO 
categories (Figure 2). In the biological process category, cellular 
processes (wounding 66.6%; defoliation 68.5%) and metabolic 
processes (wounding 61.5%; defoliation 62.5%) were the most 
dominant groups, followed by stimulus responses (wounding 
23.9%; defoliation 19.4%) and cellular component organization 
(both approximately 17%). Regarding the molecular function 
category, 65.9% of the wounding unigenes and 63.5% of the 
defoliation unigenes were assigned to binding, followed by 
catalytic activity (both approximately 49%), transcription regula- 
tory activity (4.6% for both), and structural molecule activity 
(wounding 4.0%; defoliation 7.5%). In the cellular component 
category, cell and cell part genes (wounding 72.1% and defoliation 
75.5% for both) dominated, followed by organelle (wounding 
49.9%; defoliation 52.4%) and organelle part genes (wounding 
25.5%; defoliation 26.4%). 

To reveal significandy enriched GO terms in DEGs compared 
to the reference transcriptome, GO enrichment analysis of the 
functional significance was performed. The GO terms with p< 
0.05 were considered significantly enriched among the DEGs. 
This analysis allowed us to determine the major molecular 
functions, biological processes, and cellular components that were 
represented among the DEGs (Figure S2). For the molecular 
functions, the significantly enriched GO terms of the DEGs for 
wounding and defoliation included DNA binding, serine-type 
endopeptidase inhibitor activity, oxidoreductase activity, and 
structural ribosome constituents specific to defoliation. For the 
enriched biological processes, there were 53 GO terms of DEGs 
for wounding and 40 GO terms of DEGs for defoliation, which 
represented cellular component organization and cellular compo- 
nent biogenesis (e.g., chromatin assembly, protein-DNA complex 
assembly, and nucleosome assembly), metabolic processes (fatty 
acid and oxylipin biosynthesis), lipid localization, and responses to 
stimuli (e.g. response to water and defense response to insects). For 
the cellular components, significantly enriched GO terms of DEGs 
for wounding and defoliation included chromatin, nucleosome, 
nucleus, extracellular region, and ribosome, and cell wall was 
specific to defoliation. 

KOG and KEGG Pathway Enrichment Analysis of DEGs 

A total of 832 wounding-responsive DEGs and 1,605 defolia- 
tion-responsive DEGs were classified functionally into 24 clusters 
of KOG. KOG enrichment analyses were performed by 
hypergeometric distribution testing, and the Bonferroni correction 
was applied. Significantly enriched functional clusters were 



selected based on a q-value £0. 05. The results demonstrated that 
genes associated with chromatin structure and dynamics, energy 
production and conversion, and amino acid transport and 
metabolism were over-represented in the KOG classification of 
wounding (Table 1), and chromatin structure and dynamics, 
carbohydrate transport and metabolism, translation, ribosomal 
structure and biogenesis, and the cytoskeleton were over-repre- 
sented in the KOG classification of defoliation (Table 2). 

A total of 460 wound-responsive DEGs and 881 defoliation- 
responsive DEGs were classified into the KEGG pathway. For 
wounding, the significandy enriched KEGG pathways included 
valine, leucine, and isoleucine degradation; linoleic acid metabo- 
lism; ascorbate and aldarate metabolism; starch and sucrose 
metabolism; fatty acid elongation; arginine and proline metabo- 
lism; and alpha-linolenic acid metabolism (Table 3). For defolia- 
tion, starch and sucrose metabolism; ribosome; alpha-linolenic 
acid metabolism; phenylpropanoid biosynthesis; flavonoid biosyn- 
thesis; stilbenoid, diarylheptanoid and gingerol biosynthesis; 
linoleic acid metabolism; and phenylalanine metabolism were 
over-represented in the KEGG pathway classification (Table 4). 

Enrichment of KEGG Pathways Common to Wounding 
and Defoliation 

From the above KEGG pathway enrichment analysis of the 
DEGs, we observed that linoleic acid metabolism, alpha-linolenic 
acid metabolism, and starch and sucrose metabolism were 
enriched KEGG pathways common to both wounding and 
defoliation. In the alpha-linolenic acid metabolism and linoleic 
acid metabolism pathways, we identified genes encoding lipoxy- 
genase (LOX), allene oxide synthase (AOS), allene oxide cyclase 
(AOC), 1 2-oxophytodienoic acid reductase, acyl-CoA oxidase 
(ACOX1, ACOX3), acetyl-CoA acyltransferase 1, and enoyl-CoA 
hydratase/3-hydroxyacyl-CoA dehydrogenase, all enzymes in- 
volved in JA biosynthesis. The gene expression levels were 
significandy increased for wounding and defoliation at 2 h 
(Figure 3). In starch and sucrose metabolism, genes involved in 
sucrose metabolism, including genes encoding sucrose synthase, 
fructokinase, and beta-fructofuranosidase (cell wall invertase), were 
down-regulated, while genes related to starch biosynthesis and 
metabolism and trehalose biosynthesis, such as genes encoding 
starch synthase, starch phosphorylase, beta-amylase, and trehalose 
6-phosphate synthase, were up-regulated (Figure 3). 

Defoliation-enriched KEGG Pathways 

For defoliation, the specifically enriched KEGG pathways were 
ribosome; phenylpropanoid biosynthesis; flavonoid biosynthesis; 
stilbenoid, diarylheptanoid and gingerol biosynthesis; and phenyl- 
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Figure 3. Expression levels of wounding- and defoliation- 
responsive DEGs with commonly enriched KEGG pathways, 
including linoleic acid metabolism, alpha-linolenic acid metab- 
olism, and starch and sucrose metabolism. Green to bright red: 
down-regulation to up-regulation. W, mechanical wounding; D, 
defoliation. 

doi:1 0.1 371 /journal.pone.0089495.g003 

alanine metabolism. Most ribosomal genes were down-regulated 
after defoliation at 2 h, but some were up-regulated at 24 h 
(Figure 4A). The expression levels of the genes involved in 



phenylpropanoid biosynthesis; flavonoid biosynthesis; stilbenoid, 
diarylheptanoid and gingerol biosynthesis; and phenylalanine 
metabolism were predominantly decreased (Figure 4B); these 
down-regulated genes included those encoding cinnamyl-alcohol 
dehydrogenase, peroxidase, caffeoyl-CoA O-methyltransferase, 
beta-glucosidase, phenylalanine /tyrosine ammonia-lyase, shiki- 
mate O-hydroxycinnamoyltransferase, chalcone synthase, antho- 
cyanidin reductase, and flavonoid 3',5'-hydroxylase. However, 
some genes involved in these pathways, including the genes 
encoding trans-cinnamate 4-monooxygenase, coumaroylquinate 
(coumaroylshikimate) 3'-monooxygenase, phenylalanine ammo- 
nia-lyase, and primary-amine oxidase, were up-regulated 
(Figure 4B). 

Confirmation of Read-mapped Genes by qRT-PCR 

To confirm the reliability of the Solexa/Illumina sequencing 
technology, fifteen genes were selected as described in the 
Materials and methods for quantitative RT-PCR assays. The 
expression patterns for thirteen genes were consistent between the 
qRT-PCR and the DGE analyses (Figure 5). The corresponding 
primers are listed in Table SI. The inconsistent expression 
between the qPCR and DGE analyses for the two remaining genes 
was most likely because DGE was more sensitive in detecting low- 
abundance transcripts and small changes in gene expression than 
qRT-PCR. 

Effect of Defoliation on Regrowth of Elongating Leaves 

Eight-week-old plants were defoliated to a stubble height of 
5 cm. The regrowth of the elongating leaves emerging above the 
cutting level was measured (Figure 6). The regrowth of elongating 
leaves was significantly restrained at 6 h after defoliation 
compared with the growth of the control plants, and the inhibition 
continued through the following seven days of measurement. 

Discussion 

Distinct Global Gene Transcription Changes for 
Sheepgrass under Wounding and Defoliation Treatments 

In this study, a transcriptome profiling analysis was performed 
to identify genes that are differentially expressed in wounding- and 
defoliation-treated plants. A sequencing depth of 4.8 to 6.4 million 
reads per library was reached (Table S2). The results of the 
saturation sequencing analyzed in the seven libraries indicated that 
the sequencing depth was sufficient to cover the transcriptome 
(Figure SI). We identified 1,836 and 3,238 wounding-responsive 
and defoliation-responsive DEGs, respectively, which suggested 
that defoliation led to nearly twice as many DEGs as wounding. 
Among these DEGs, 954 were uniquely expressed under 
defoliation and thus represent candidate genes that merit further 
study. Further analysis indicated that the defoliation-responsive 
gene number was highest at 2 h, during which the number of 
down-regulated defoliation-responsive DEGs (1,526) were over 
four times the number of down-regulated wounding-responsive 
DEGs (350), suggesting that plants could rapidly respond to the 
enormous damage of defoliation by regulating global gene 
expression. 

Wounding and Defoliation Commonly Activate Systemic 
JA Synthesis 

JA is derived from the unsaturated fatty acid a-linolenic acid 
(18:3), an abundant octadecanoid in the cellular membranes of 
higher plants. Lipoxygenase (LOX) is the initial step in JA 
biosynthesis and catalyzes the oxygenation of oc-linolenic acid, 
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Figure 4. Expression levels of defoliation-responsive DEGs specifically enriched in KEGG pathways. (A) Ribosomal DEGs. (B) DEGs 
involved in phenylpropanoid biosynthesis; flavonoid biosynthesis; stilbenoid, diarylheptanoid and gingerol biosynthesis; and phenylalanine 
metabolism. Green to bright red: down-regulation to up-regulation. D, defoliation. 
doi:1 0.1 371 /journal.pone.0089495.g004 



leading to the formation of hydroperoxy octadecadienoic acids, 
which are then converted into 12-oxo phytodienoic acid (12- 
OPDA) in the plastids via enzymatic reactions catalyzed first by 



allene oxide synthase (AOS) and subsequently by allene oxide 
cyclase (AOC). OPDA is then converted to JA by 12-oxo- 
phytodienoic acid reductase (OPR) followed by three cycles of (3- 
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Figure 5. Quantitative PCR validations of read-mapped genes. 

doi:1 0.1 371 /journal.pone.0089495.g005 



oxidation in the peroxisome [26,27]. In this study, we observed 
that the DEGs involved in JA biosynthesis after wounding and 
defoliation were significandy enriched according to the GO terms 



and KEGG pathway enrichment analysis (Figure S2, Table 3, 4). 
The transcription levels of key genes in JA biosynthesis, including 
genes encoding lipoxygenase (LOX), allene oxide synthase (AOS), 



9 




Time after defoliation (h) 

Figure 6. Impact of defoliation on regrowth of elongating leaves above the cutting level. After 8 weeks of growth, sheepgrass plants 
were defoliated to a stubble height of 5 cm, and regrowth of elongating leaves was measured above the cutting level after defoliation at different 
time points. 

doi:1 0.1 371 /journal.pone.0089495.g006 
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allene oxide cyclase (AOC), 12-oxophytodienoic acid reductase 
(OPR), acyl-CoA oxidase (ACOX1, ACOX3), acetyl-CoA acyl- 
transferase 1, and enoyl-CoA hydratase/3-hydroxyacyl-CoA 
dehydrogenase (Figure 3), were significandy increased after 
wounding and defoliation. Our results are consistent with the 
previous reports, demonstrating that wound-induced responses are 
triggered by the de novo synthesis of the plant hormone JA when 
plant tissues are injured by various means, including herbivores, 
pathogens, and mechanical stress [9,28-29]. 

Defoliation Specifically Reduces the Expression Levels of 
Ribosomal Protein Genes, Cell Division or Cell Expansion- 
related Genes, and Lignin Biosynthesis Genes and may 
Negatively Affect Plant Growth 

Enrichment analyses of GO terms, KOG, and the KEGG 
pathway suggested that the ribosome was specifically over- 
represented after defoliation. The ribosome is an essential 
ribonucleoprotein complex that is engaged in translation and is 
indispensable for growth [30]. In this study, many ribosomal 
protein genes, including the large 60S and small 40S subunits of 
the cytosolic ribosome, were identified, and their expression levels 
mainly decreased at 2 h after defoliation (Figure 4A). In both 
plants and animals, the loss of ribosomal proteins leads to reduced 
growth, which most likely correlated with reduced ribosome 
production and lower rates of protein synthesis [31]. For example, 
in Arabidopsis thaliana, several RP loss-of-function mutations that 
affect cell division or cell expansion and consequendy result in a 
deformed leaf size and shape have been identified, indicating cell- 
or development-specific roles of RPs during leaf growth [32]. 
Notably, in this study, the expression levels of numerous genes 
involved in cell division or cell expansion were down-regulated 
after defoliation, and these genes include those encoding cyclins, 
cyclin-dependent kinases, and expansion (Table S4). 

Lignin is a phenolic heteropolymer of the secondary cell walls 
that plays a major role in plant development and defense against 
pathogens. The monolignols represent the main component of 
lignin, and their synthesis involves many intermediates and 
enzymes. In this study, we identified several key genes involved 
in lignin biosynthesis that were down-regulated, including genes 
encoding phenylalanine ammonia-lyase (PAL), caffeoyl-CoA O- 
methyltransferase (CCOMT), shikimate O-hydroxycinnamoyl- 
transferase (HCT), and cinnamyl alcohol dehydrogenase (CAD) 
(Figure 4B). The down-regulation of these genes has been 
demonstrated to reduce lignin biosynthesis and cell-wall thickness, 
alter xylem organization, and retard growth [33-36]. 

Notably, the regrowth of elongating leaves was significandy 
retarded after defoliation at 6 h through the following 7 days of 
measurement (Figure 6), suggesting that the above gene expression 
patterns correlated consistently with the phenotype. Castrillon- 
Arbelaez et al. demonstrated that although defoliation led to a 
rapid and transient reduction of non-structural carbohydrates 
(NSC) in grain amaranth, only a few changes in gene expression 
and enzyme activity could be associated with the NSC changes in 
an analysis of 25 genes. They also revealed that the rapid 
mobilization of foliar starch reserves followed by an efficient 
recovery of all NSC reserves after defoliation did not alter plant 
growth or reproductive fitness [17]. Vargas-Ortiz et al. demon- 
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